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Abstract 

A general turbulent constitutive relation (Shih and Lumley, 1993) is directly applied 
to propose a new Reynolds stress algebraic equation model. In the development of this 
model, the constraints based on rapid distortion theory and realizability (i.e. the positivity 
of the normal Reynolds stresses and the Schwarz’ inequality between turbulent velocity 
correlations) are imposed. Model coefficients are calibrated using well-studied basic flows 
such as homogenous shear flow and the surface flow in the inertial sublayer. The performance 
of this model is then tested in complex turbulent flows including the separated flow over a 
backward- facing step and the flow in a confined jet. The calculation results are encouraging 
and point to the success of the present model in modeling turbulent flows with complex 
geometries. 

1. Introduction 

The present study concentrates on complex turbulent shear flows which are of great 
interest in propulsion systems. The particular flows presented in this paper are for the 
backward-facing step and the confined jet, both of which have complex structures. For 
example, a confined jet combines several types of flow structure and flow phenomena such 
as a shear layer, jet, recirculation, separation and reattachment. Accurate prediction of 
these flows is of great importance in all the key elements of engine design. 

The turbulence model developed in this study is a Reynolds stress algebraic equation 
model which is based on a turbulent constitutive relation (Shih and Lumley, 1993), a result of 
rapid distortion theory (Reynolds, 1987) and the turbulent realizability principle (Schumann 
1977, L uml ey, 1978). The constitutive relation is obtained using the invariance theory in 
continuum mech ani cs. For flows including a passive scalar, this theory leads to a general 
constitutive relation for the Reynolds stress tensor ujuj in terms of the mean deformation 
rate tensor U it j and the turbulent velocity and length scales characterized by the turbulent 
kinetic energy k and its dissipation rate e. Pope (1975) applied a similar constitutive 
relation to Rodi’s algebraic Reynolds stress formulation (Rodi, 1972) in conjunction with 
the LRR second order closure model (Launder et al., 1975) and obtained an explicit algebraic 
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expression for the Reynolds stresses for a two-dimensional mean flow field. Taulbee (1992, 
1994) was able to extend this method to a general three-dimensional flow. Gatski and 
Speziale (1992) also applied this method in their algebraic Reynolds stress model. We note 
that in Rodi’s algebraic Reynolds stress formulation, some assumptions, such as the constant 
anisotropy of Reynolds stresses and the neglect of turbulent transport of second moments, 
are in general not valid for most turbulent shear flows. These assumptions may bring large 
errors to turbulence modeling. In addition, the deficiency of existing second order closure 
models would also add extra errors to this type of model. In this study, Rodi’s formulation 
is not considered. As an alternative, we directly impose the constraints based on rapid 
distortion (rotation) theory and realizability on the constitutive relation for the Reynolds 
stresses. As a result, a realizable algebraic expression for the Reynolds stresses in terms of 
the mean velocity gradient and the characteristic scales of turbulence is obtained for general 
three-dimensional turbulent flows. For turbulent scales, the standard k-e model transport 
equations are used in this study. Some model constants are calibrated using a well-studied 
homogeneous shear flow and a surface flow in the inertial sublayer and then tested in other 
complex flows. The model validation is made on the basis of applications to the rotational 
homogeneous shear flows simulated by Bardina et al. (1983), the two backward-facing step 
flows experimentally studied by Driver and Seegmiller (1985) and Kim et al. (1978) and the 
five cases of confined jets studied by Barchilon and Curtet (1964). 

The calculations for complex flows are performed with a conservative finite volume 
method (Zhu, 1991b). Grid independent and low numerical diffusion solutions are obtained 
by using differencing schemes of second-order accuracy on sufficiently fine grids. For wall- 
bounded flows, the standard wall function approach (Launder and Spalding, 1974) is used 
for wall boundary conditions. The results are compared in detail with the experimental 
data for both mean and turbulent quantities. The calculations using the standard k-e eddy 
viscosity model are also carried out for the purpose of comparison. The comparison shows 
that the present realizable Reynolds stress algebraic equation model significantly improves 
the predictive capability of k-e equation based models, expecially for flows involving massive 
separations or strong shear layers. In these situations, the standard eddy viscosity model 
overpredicts the eddy viscosity and, hence, fails to accurately predict shear stress, adverse 
pressure gradient, separation, reattachment, recirculation, etc. We find that the success of 
the present model in modeling complex flows is largely due to its effective eddy viscosity 
formulation which accounts for the effect of the mean deformation rate. According to 
the present model, the effective eddy viscosity will be significantly reduced by the mean 
deformation rate and maintained at a correct level to mimic the complex flow structures. 
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2. Turbulence model 


2.1 Constitutive relation. Constitutive relations for the Reynolds stresses were 
derived by several researchers (Pope, 1975, Yoshizawa, 1984 and Rubinstein and Barton, 
1990). Shih and Lumley (1993) used the invariant theory in continuum mechanics and 
the generalized Cayley-Hamilton formulations (Rivlin, 1955) to derive a more (perhaps the 
most) general constitutive relation for the Reynolds stresses under the assumption that the 
Reynolds stresses are dependent only on the mean velocity gradients and the characteristic 
scales of turbulence characterized by the turbulent kinetic energy k and its dissipation rate 
e. This relation is 


UiUj = -zkSij + 2a 2 {Ui,j + Uj,i ~ + 204 (U it j + Uj i glli^jj) 


K 3 


+ 2 a^{U itk Uj,k - ilMii) + 2a^{U kti U kyj - ^n 2 5 0 ) 

+ 2 a s ^-(Ui, k Ul k + Ul k U jtk - ^n 3 5 0 ) + 2 a x ^{U k ,iUli + ~ 

+ 2 a 12 ^-(Ul k Ul k - ^Il^i) + 2a 13 ^-(Ul i Ul J - |ll4«y) 

+ 2a 14 ^ r (u i , k u l<k ul j + u hk u Uk ul - |n 5 5ii) 

+ 2 a 16 ¥-(u i , k ul k u? tj + u j>k ul k uf ti - |n.««) 

+ 2 a lt ¥-(U itk U ltk Ul m Ul m + U jlk Ui, k Uf tn JJf tm - iSij) 


(1) 


where 


Hi UtykUkytl 

n 4 = ultVU* 


n 2 = Ui, k Ui, k , n 3 = u itk ul k , 
n 5 = Ui tk Ui <k u? ti , n 6 = Ui, k uf tk ul i} 


n 7 = Ui,kU ltk ul n ul r 


( 2 ) 


Eq.(l) contains 11 undetermined coefficients which are, in general, scalar functions of various 
invariants of the tensors in question, for example, SijSij (strain rate) and ClijClij (rotation 
rate) which are (H 2 +Hi)/2 and (R 2 -Hi)/2 respectively. The detailed forms of these scalar 
functions must be determined by other model constraints such as rapid distortion theory, 
r ealiz ability, and appropriate experimental data. 

It is noticed that the standard k-t eddy viscosity model corresponds to the first two 
terms on the right hand side of Eq.(l). Both the two-scale DIA approach (Yoshizawa, 1984) 
and the RNG method (Rubinstein and Barton, 1990) also provided a similar relation which 
is the first five terms on the right hand side of Eq.(l). 
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Iii this study, for the purpose of engineering application we truncate Eq.(l) to its 
quadratic tensorial form. The necessity of using higher order non-linear terms will be left 
for future study. To distinguish between the strain and rotation, we define 


where 


stj = Sij - -SuSij , % = 

c(2*) __ q2 ^ q2 r o< 2 *> — O 2 ^0 2 A • 
b ij — *ij ~ > ll ij — “ ^ L kk d tJ 


Sij = + Uj,,), fly = i (Uij - Uj,i ) 


= SikSkj tiij = Clik^lkj 

From the above definitions, we have the following relations: 

Sh = 0, 4 2,) = 0, s4 2-) = o 

For later use, we further define 


C* Q* C* c(2*) c(2*) c(2*) 

pp. ^ij^jk^ki ti/-( 2*) _ gjj ^ 

(S*) 3 (S'C 2*))3 

tf* = , p ( 2,) = if* + ng ,) ng* ) 


(3.1) 

(3.2) 


(3.3) 


(3.4) 


Using Eqs.(3.1-3.4), the truncated equation (1) can be written as 


mj = \kfij - c»-2s; + c 1 ^2(4 2,) + «<?•>) 

+ - < -> - + nisij) 

+ - niT* + (4) 


2.2 Rapid distortion constraint. Reynolds (1987) and Mansour et al. (1991) 
studied the effect of rapid rotation on turbulence using rapid distortion theory (RDT). 
It was shown that there is no effect of the rapid mean rotation on the isotropic turbulence. 
This result provides a constraint for Eq.(4). For rotating flows with Sij = 0, Eq.(4) becomes 



UiUj 

~2k~ 



k? 


= sX’-'W - c 2 - C 3 ) 
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From the result of RDT for the isotropic turbulence, bij should remain zero under rapid 
mean rotation and, therefore, we must require 2Ci = C 2 + C 3 - As a result, Eq.(4) becomes 


Ui Uj = - Cp — ZSij 


+ c^(2s^ ) -s- k n- tl + n-„s; j ) 

+ C ^( 2S <?' + s * n h - 


(S) 


2.3 Re aliz ability. Realizability (Schumann, 1977, Lumley,1978), defined as the re- 
quirement of the non-negativity of turbulent normal stresses and Schwarz’ inequality be- 
tween any fluctuating quantities, is a basic physical and mathematical principle that the 
solution of any turbulence model equation should obey. It also represents the minimal 
requirement to prevent a turbulence model from producing unphysical results. In the fol- 
lowing, this principle will be applied to the relation of Eq.(5) to obtain constraints on its 
coefficients C 2 and C 2 . The same procedure together with the RDT constraint can be 
also applied to the full equation (1). 

Turbulence models often produce unphysical results under some extreme situations. For 
example, under a rapid mean strain the turbulent energy component in the strain direction 
will be rapidly reduced and a non-realizable model often drives that energy component to a 
negative value and under a high mean shear the turbulent shear stress will rapidly increase 
and a non-r ealiz able model often overpredicts this increase such that the Schwarz’ inequality 
will be violated. The commonly used fc-e eddy viscosity model with a constant = 0.09: 


UiUj 



( 6 ) 


is one such unrealizable model. In this model, the energy component u\ will become negative 
when S^k/e > 1/0.27 and the correlation coefficient between u 2 and u 2 will exceed unity 
when - > 1/0.27 for a pure mean shear flow (which has only one non-zero component 

s; 2 ). 

To make eddy viscosity model Eq.(6) realizable, the coefficient cannot be a constant. 
It must vary with the mean flow deformation rate. To determine its appropriate formulation, 
we may use the following realizability constraints: 

^>0 (a = 1,2,3) 

§£<1 (a = 1,2,3; (3 = 1,2,3) 

U l 
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(7.1) 

(7.2) 



Reynolds (1987) used the constraint of Eq.(7.1) to formulate the coefficient which en- 
sures positive normal stresses. Shih et al. (1993) also imposed Eq.(7.1) on their Reynolds 
stress algebraic equation model. Here, we will follow the method of Reynolds and use both 
Eqs.(7.1) and (7.2) to determine the coefficients in Eq.(5). 

In the principal axes of S*j (note, 5,* = 0), we may write: 



0 

1+0 

2 

0 


The invariant S* and W* defined in Eq.(3.4) can 


0 

0 


s : 


ii 


_i=« i 
2 ' 

be calculated as 


s- = |s;,| 


*±* w = 


( 8 ) 


In addition, noting that in the principal axes of S*j the off-diagonal terms of 5^ 2 *^ are also 
/2*\ 

zero and that S). ’ = 0, we may write 


q(2*) _ 
^ij 


The invariants S^ 2 *^ and W^ 2 *^ are 


' 1 
0 

,0 


0 

l±fc 


0 


0 

0 

1-J> 


? ( 2 *) 

'11 




( 9 ) 


According to Eq.(5), the energy component u 2 in the principal axes of Stj is 

ui=h- 25^ + (C 2 + C' 3 )^25S?* ) 

Now let us consider the contraction case in which S > 0, 5”^*^ > 0 and u\ will decrease 
due to the contraction strain. Using Eqs.(8) and (9), we obtain 

s* - c '“T 2 S '\f^ + ( ° 2 + ft £ 2S<!-) \/rS 

Now, applying the constraint Eq.(7.1) and allowing the component u\ — * 0 but remains 
positive as S* — * oo and S^ 2 *^ — » oo. To satisfy this constraint, we may let 

A 


cv = 


C 2 + C*3 = 


+ \[ 


1L -U*k 

3 +a* U e 
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and 


A- B = 1 


Following Reynolds (1987), let 


K 


18 


3 + a 2 


A? 


^ = t 


18 


3 + 6 2 


( 10 ) 


Using Eqs.(8) and (9), A* and A? m) can be determined by the foUowing equations: 

(4;)J_?i;_9r=0 (11.1) 

(^(2*))3 _ £ A (2*) _ 9W t(2.) = 0 (11.2) 

It can be shown that the positive root of the above equations can be obtained when the 
values of W* and W < 2 *> are between -l/y/6 and l/y/6 which correspond to axisymmetric 
expansion and axisymmetric contraction respectively. The appropriate roots are 

A* = V6 cos 4>, <f> = \arc cos (V^W*) (H-3) 

o 

A^ 2 *^ = \/6cos<£, <f> = — arc cos (H*4) 

o 

Eqs.(11.3) and (11.4) show that the values of A* and A? m) are between y/6/2 and y/6. The 
model coefficients can be now written as 

c = A 
M A 0 + A* t U*± 

B 

° 2 + Ci ~ A> + 

The further determination of A, B, Ao and A\ should be carried out by using the constraint 
of Eq.(7.2) and the experimental data from well-studied turbulent flows such as homogeneous 
shear flows and channel flows. 

Here, we try to propose a simple as possible but workable model (which contains the 
property of anisotropy) for engineering application and leave the more complete model form 
of Eqs.(5), (12) and (13) for future study. To do that, we choose .4 = 1, then B must be 
equal to zero and C% = — C 2 . As a result, Eq.(5) becomes 

u-uj = \kSij - <?„— 2S£ + 2 C 2 ^{-St k nij + WkS'kj) (14) 

3 £ ^ 


( 12 ) 

(13) 
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It is obvious that this model satisfies the constraint Eq.(7.1). To apply the constraint 
Eq.(7.2), we use a pure shear flow with only one non-zero component Ui }2 (i.e., 5* 2 = flj 2 > 
0) which can be considered as the most extreme case for satisfying Schwarz’ inequality. For 
this flow, the relevant Reynolds stresses are 


— C^ 2jl>i 2 

^ = ^ + 4C 2 ^s: 2 nt 2 (15) 

o £* 

^ = §*-4C 2 ^s; 2 n; 2 

Now using the constraint of Eq.(7.2), we may find a formulation for C 2 : 


C 2 


yA - 9Cl{S±y 


(16.1) 


where 


C u = 


A 0 + A 


(16.2) 


The model represented by Eqs.(14), (16.1) and (16.2) is quite simple but has several ad- 
vantages compared to the standard k-e eddy viscosity model of Eq.(6). First, the present 
model is fully realizable. It will not produce negative energy components and will not violate 
the Schwarz’ inequality between turbulent velocities. Second, the effective eddy viscosity, 
defined as u^ap/2 5*^, is anisotropic as it should be. Finally, the present model contains 
the effect of mean rotation on Reynolds stresses with a proper behavior that matches the 
RDT result: the mean rotation will not affect the isotropic turbulence. 

There are still two model constants, Aq and Co, that need to be determined. We may 
use Eqs.(15) for the homogeneous shear flow or the surface flow in the inertial sublayer. 
According to these flows, Aq and Co are chosen as 


A 0 = 6.5, Co = 1.0 


(17) 


With the values of Ao and Co in Eq.(17), the model of Eq.(14) gives 6 i 2 = —0.156, = 

— b 22 = 0.123 for Tavoularis and Corrsin’s (1981) homogeneous shear flow at U\ >2 k/e = 6.08 
and gives b\ 2 — —0.122, 5u = —b 22 = 0.14 for the direct numerical simulation of channel flow 
(Kim, 1990) in the inertial sublayer at U\ f2 k/e = 3.3. These results show that the present 
model gives reasonable anisotropy of Reynolds stresses for both the homogeneous shear flow 
and the boundary layer flow compared to the standard k-e eddy viscosity model which gives 
5n = b 22 = 0 for both the flows and gives b\ 2 = —0.273 for the homogeneous shear flow and 
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l 12 — -0.149 for the boundary layer flow. Detailed comparisons with the experimental and 
DNS data are shown in Table 1 for the homogeneous shear flow of Tavoularis and Corrsin 
(1981) and in Table 2 and Figure 1 for the channel flow of Kim (1990). 


Table 1. Anisotropy in the homogeneous shear flow 



experiment 

standard 

present 

i>12 

-0.142 

-0.273 

-0.156 

fell 

0.202 

0 . 

0.123 

fe 22 

-0.145 

- - 

0 . 

-0.123 


Table 2. Anisotropy in the channel flow 



DNS data 

standard 

present 

fel2 

-0.145 

-0.149 

-0.122 

feu 

0.175 

0 . 

0.14 

fe 22 

-0.145 

0 . 

-0.14 


2.4 Model equations. Here we summarize the equations and the models which 
will be used for applications in the next section. For incompressible flows, the mean flows 
are governed by the following equations 


Ui,i = 0 

Vi,, + (VjUi - + vTi uj)j = 

where the Reynolds stresses will be modeled by Eq.(14): 

U{Uj = -kSij — 0 ^ — 2 S*j + 2C?2— (— + flSfcSy 

and (7„, C'j are determined by Eq.(16): 


(18) 

(19) 


c *-3TMs=r C2 = 


Vl - 9 CJ (^) 2 


C 0 +6^^ k 


where 


A 0 = 6.5, Co = 1.0 

Two quantities in Eq.(14), the turbulent kinetic energy k and its dissipation rate e, remain 
to be determined. At the present time, we use the standard k-e model equations which are 

k,t + Ujkj = \{v + — )fc,j],j — uvjUij — £ (20) 

£,t + Uj£,j = [{u + - C'l^UiUjUij - C'l— (21) 
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( 22 ) 


(23) 

3. Applications 

3.1 Rotating homogeneous shear flow. The present model is able to mimic the 
effect of the mean rotation rate on the turbulence. A test case is the rotating homogeneous 
shear flow which was studied by Bardina et a J. (1983) using the large eddy simulation 
(LES) method. The effect of solid body rotation or the rotation of the reference frame on 
the turbulence must be appropriately incorporated in Eq.(14) through the terms containing 
In addition, the coefficients and C 2 should be also modified by the rotation rate of 
the reference frame, u>i (angular velocity). Particularly, the U* in Eq.(16.2) is modified by 

V = yjSfjSl, (24.1) 

where 

= f l*j - 2eijkWk 

ij = Clij Gijk^lc 

where Clij is the mean rotation rate viewed in the rotating reference frame. Figure 2 is 
the configuration of the flow being tested where Cl — u 3 and = Sf 2 = | dU/dy — 
S/2. Figures 3(a)-3(c) show the evolution of the turbulent kinetic energy Jc/ko with the 
nondimensional time, 5f, at the rotation rates of Cl/S = 0, 0.5 and —0.5, respectively, 
where ko is the initial turbulent kinetic energy, S is the mean strain rate and Cl is the 
angular velocity of the reference frame. The calculations were performed with a fourth 
order Runge-Kutta scheme. The initial condition corresponding to the isotropic turbulence 
used in LES with £0 f Sk 0 = 0.296 was adopted for all the three cases. The results from both 
the present model and the standard k-e model (hereafter referred to as SKE) are compared 
with LES results in figures 3(a)-3(c). These figures show the ability of the present model 
to simulate the effect of the large rotation rate on turbulence. Note that the SKE model 
gives the same results as for the no rotation case because it cannot account for the effect of 
rotation on the evolution of turbulence. 


(24.2) 


where 

u t = Cp — 

£ 

The coefficients C e i, C e 2 ,<Tk and <j e assume their standard values: 

<7 e i = 1.44, C e 2 = 1.92, o k = 1, <r e = 1*3 
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3.2 Backward-facing step flows. 

Numerical procedure. For computational convenience, the non-dimensional form 
of the governing equations is solved, in which 


< Xi > 



< Ui >= 



<p>= 



< k > = 



< e>= 


£L re f 



< Vt >= 


Vt 


Uref L re f 


(25) 


where < > refers to a non-dimensional quantity, and L r cfi and U re f are the reference length 
and velocity, respectively. Accordingly, the flow Reynolds number is defined by 


Re = LrefU ^ l (26) 

V 

Hereafter, all the quantities will be of the non-dimensional form so that < > will be dropped 
for simplicity. 

In the steady-state and two dimensional cases (xi = x,®2 = y)> the transport equations 
(19), (20) and (21) can be written in the following general form 

[U4> ~^Te + + [V<f> + 7+ )4> ’ v] ' V = S * (27) 

where <f> stands for the dependent variables: U, V , k and e. S# is the source term for each 
corresponding equation. 

The n umer ical method used to solve the system of equations (27) is a finite-volume 
procedure. It uses a non-staggered grid with all the dependent variables being stored at the 
geometric center of each control volume (Figure 4). The momentum interpolation procedure 
of Rhie and Chow (1983) is used to avoid spurious oscillations usually associated with 
the non-staggered grid, and the pressure-velocity coupling is handled with the SIMPLEC 
algorithm (Van Doormal and Raithby, 1984). To ensure both accuracy and stability of 
the numerical solution, the convection terms are approximated by a second-order accurate 
and bounded differencing scheme (Zhu, 1991a), and all the other terms by the conventional 
central differencing scheme. As a result, the discretized counterpart of equation (27) can be 
cast into the following linearized form 


4>c ^2, Ai — A[4n + Sc 

i i 


(28) 


where the coefficients Ai (l = W,E,S,N), which relate the principal unknown 4>c to its 
neighbours <f>i (Figure 4), result from the discretization of the left-hand side terms of equation 
(27). The convection scheme used ensures that Ai > 0 so that the resulting coefficient matrix 
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is always diagonally dominant. The strongly implicit procedure of Stone (1968) is used to 
solve the system of algebraic equations. The iterative solution process is considered to be 
converged when the maximum normalized residue of all the dependent variables is less than 
10 -4 . The details of the present numerical procedure are given in Rodi et a J. (1989) and 
Zhu (1991b). 

Numerical results. Application is made to the two backward-facing step flows ex- 
perimentally studied by Kim, Kline and Johnston (1978) and Driver and Seegmiller (1985), 
from here on referred to as KKJ- and DS-cases, respectively. Figure 5 shows the flow con- 
figuration and the Cartesian coordinate system. Table 3 gives the flow parameters for both 
cases; here the experimental reference free-stream velocity U re f and step height H t are taken 
as the reference quantities for non- dimensionalization. 


Table 3. Flow parameters 


Case 

Re 

S 

L s 

L e 

H, 

H d 

Uref 

DS 

37423 

1.5 

10 

40 

1 

8 

1 

KKJ 

44737 

0.6 

10 

40 

1 

2 

1 


Three types of boundaries are present, i.e. inlet, outlet and solid wall. At the inlet, 
the experimental data sue available for the streamwise mean velocity U and the turbulent 
normal stresses uu and vv. k is calculated from these uu and vv with the assumption that 


ww = 


-{uu + vv) 


(29) 


and £ by 

(J 3 / 4 1 . 3/2 

e = , L = min(0.41Ay, 0.0855) (30) 

JLf 

where Ay is the distance from the wall and S is the boundary-layer thickness given in Table 
3. At the outlet, the streamwise derivatives of the flow variables are set to zero. Influences 
of both inlet and outlet conditions on the solution are examined by changing the locations 
L a and L e , and it has been found that in both cases, the distances given in Table 3 are 
already sufficiently far away from the region of interest. In the earlier stage of this work, 
we tested several low Reynolds number k-e models including those of Chien (1982), Lam 
and Bremhorst (1981), Launder and Sharma (1974), Shih and Lumley (1992), and Yang 
and Shih (1992), but none of them was found to be able to yield satisfactory solutions for 
the skin friction along the bottom wall. Similar findin g s were also reported in Awa et a J. 
(1990), Shuen (1992) and So and Lai (1988). Therefore in this work, we use the standard 
wall function approach (Launder and Spalding, 1974) to bridge the viscous sublayer near 
the wall. 

Two sets of non-uniform computational grids are used to examine the grid dependence 
of the solution; they contain 110x52 (coarse) and 199x91 (fine) points for the KKJ-case and 
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106x56 (coarse) and 201x109 (fine) points for the DS-case. Figures 6(a) and 6(b) show the 
friction coefficient Cf at the bottom wall calculated with the SKE model and the present 
model; also included in figure 6(a) are the experimental data for the DS-case, but no such 
data are available for the KKJ-case. It can be seen that the grid refinement does produce 
some differences for the results of the present model, more noticeable in the KKJ-case, and 
this is also the case for the SKE results. This indicates that the solutions obtained on the 
coarse grids are not sufficiently close to the grid-independent stage. Recently, Thangam 
and Hur (1991) have conducted a highly-resolved calculation for the KKJ-case. They have 
found that quadrupling a 166x73 grid leads to only a minimal improvement. Therefore, 
the present results on the fine grids can be considered as grid-independent. For the DS- 
case, the fine grid computations with the SKE model and present model required 703 and 
691 iterations, and took approximately 7.1 and 9 minutes of CPU time on the Cray YMP 
computer. In the following, only the fine grid results are presented. 

The wall friction coefficient Cf is a parameter that is very sensitive to the near- wall 
turbulence modeling. It is Cf that the various low Reynolds number k-e models tested 
predict much worse than those using wall functions. However, the influence of the near-wall 
turbulence modeling is mainly restricted to the near-wall regions. It is seen from figure 6(a) 
that both the SKE model and the present model largely underpredict the negative peak of 
Cf , pointing to limited accuracy of the wall function approach in the recirculation region. 

The computed an d measured reattachment points are compared in Table 4. They are 
determined in the calculation from the point where Cf goes to zero. The reattachment 
point is a critical parameter which has often been used to assess the overall performance 
of turbulence models as well as numerical procedures. Table 4 clearly demonstrates the 
significant improvement obtained with the present model. It is important to mention that 
this im provement is mainly due to the behavior of C p in the present model, and that the 
anisotropic behavior of the turbulent stresses only makes a marginal contribution to it. 


Table 4. Comparison of the reattachment points 


Case 

measurement 

SKE 

PRESENT 

DS 

6.1 

4.99 1 

5.80 

KKJ 

7± 0.5 

6.35 

7.27 


Figures 7(a) and 7(b) show the comparison of computed and measured static pressure 
coefficients C p along the bottom wall. In both cases, the SKE model is seen to predict 
premature pressure rises which is consistent with its underprediction of the reattachment 
lengths. 

The streamwise mean velocity U profiles are shown in figures 8(a) and 8(b) at four 
different cross-sections. Here, the differences between the results of the SKE model and 
present model are not substantial, as compared to other flow variables. However, the present 
model shows somewhat slower recovery in the vicinity of the reattachment point. We notice 
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that such a slow recovery also exists in the Reynolds stress model prediction by Obi et a 1. 
(1989). Further downstream, say at ®=20 in figure 8(a), the results of the two models nearly 
coincide with each other. 

Finally, the comparisons of predicted and measured turbulent stresses u 2 , v 2 and uv 
are shown in figures 9 and 10 at various ar-locations. In the KKJ-case, no experimental data 
for the turbulent stresses are available in the recirculation region, and the reattachment 
point was found in the experiment to move forward and backward continuously around 
seven step heights downstream of the step, leaving an uncertainty of ±0.5 step height for 
the reattachment length. This also points to some uncertainty in the measured turbulent 
quantities in the recovery region. On the other hand, the experimental data in the DS-case 
should be considered more reliable because of the smaller uncertainty of the reattachment 
location, indicating a smaller unsteadiness of the flow. The SKE model gives unrealistic 
results about normal Reynolds stresses: v 2 > u 2 at all the locations. In contrast, the 
present model gives at least qualitatively correct results due to the non-linear terms in 
Eq.(14) which increase u 2 while decreasing v 2 , leading to an overall improvement in both 
u 2 and v 2 results. 

3.3 Confined Jets. The general features of confined jets measured by Barchilon and 
Curtet (1964) are sketched in figure 11. At the entrance, two uniform flows, a jet of larger 
velocity and an ambient stream of smaller velocity, are discharged into a cylindrical duct of 
diameter D 0 . The inlet flow conditions can be characterized by the Craya-Curtet number 
C t . The experiment shows that recirculation occurs when C t <0.96. For a given geometry, 
recirculation as well as adverse pressure gradients can be intensified by reducing the value of 
C t at the entrance. Five cases of Ct were studied, ranging from no to strong recirculation. 

The predicted axial mean velocity profiles at two Ct numbers axe shown and compared 
with the experimental data in figure 12, where R and Um are the radius of the cylinder and 
the sectional mean velocity, respectively. Both models are seen to predict very well the up- 
stream evolution of the flow. As for the downstream development, the results of the present 
model remain in good agreement with experiments while the SKE model underpredicts the 
centerline velocity decay at all Ct numbers. 

The variation of the pressure coefficient C p along the duct wall is shown in figure 13. 
The pressure distribution is governed by the jet entrainment as well as the contraction and 
expansion of the flow caused by the recirculation bubble. The decrease in the ambient 
velocity induced by the entrainment gives rise to an adverse pressure gradient, while the 
contraction of streamlines produces the opposite effect. These two mechanisms interact 
more intensely with each other as Ct decreases and cause the pressure to vary little in the 
region upstream of the center of the recirculation bubble. However, in the downstream part 
of the recirculation bubble, the deceleration of the flow sets up an adverse pressure gradient, 
the slope of which becomes steeper as Ct decreases. Therefore, the ability to capture the 
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location of the recirculation center will have a direct impact on the prediction of the pressure. 
Regarding the comparison between predictions and experiments, it is seen that although 
both models predict the same total pressure rises which are in excellent agreement with the 
measurements , the present model captures the pressure distribution much better than does 
the SKE model for all the C t values. 

4. Conclusion 

A new Reynolds stress algebraic equation model has been developed using a truncated 
constitutive relation. The development of the model is based on the constraints from rapid 
distortion (rotation) theory and realizability. Therefore, the present model shows the proper 
lack of a rotation effect on the isotropic turbulence and is fully realizable, i.e., it will not pro- 
duce unphysical Reynolds stresses for the mean flow field. The model is calibrated by using 
basic flows (homogeneous shear and channel flows), and then is applied to complex flows. 
The calculations have been compared with available experimental data. The comparisons 
show that the present model does provide significant improvement over the standard k-e 
eddy viscosity model and that the present model is as robust and economical as well. This 
indicates that the present model has good potential to be a practical tool in engineering 
applications. 
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Figure 4 Typical control volume centered at C and related nodes 



Figure 5 Backward-facing step geometry 
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Figure 7 Pressure coefficient 
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Figure 8 Streanrwise mean velocity U-profiles. 
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Figure 9 Turbulent stress profiles in DS-case. 
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Figure 10 Turbulent stress profiles in KKJ-case. 
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Figure 11 Flow configuration for confined jet 
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Figure 12 Axial mean velocity profiles, 
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